SNP predition ¶
Author: Bartosz Lewandowski¶
------------------------------------------------------------------------------------------------------------¶
Workflow¶
1. Understanding our data and quick review
2. Preprocessing
3. Building first neural network
4. Checking results
5. Plots for our first neural network
6. Bulding second neural network
7. Checking results
8. Plots for second NN
9. Last NN using StratifiedKfold
10. Results
11. Summarizing results and choosing best neural network Import of necessary function¶
In [1]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sn
import pandas_profiling
import tensorflow as tf
from sklearn.model_selection import train_test_split
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from keras.models import Sequential
from keras.layers import Dense
from keras.callbacks import EarlyStopping
from keras.callbacks import ModelCheckpoint
from sklearn.metrics import balanced_accuracy_score
from sklearn.metrics import matthews_corrcoef
from sklearn.metrics import roc_auc_score
from sklearn.impute import KNNImputer
from pandas_profiling import ProfileReport
from keras import backend as K
from keras.layers import Dropout
from sklearn.metrics import roc_curve,auc
1. Data overview¶
Reading data¶
In [3]:
data = pd.read_csv('data/logreg.txt', sep=';', dtype={'genotype':'float32'})
In [3]:
data.head()
Out[3]:
We are checking how many variables and observations we have. We can see that there are 12 variables, 1 dependant and 11 independant. The dependant variable is genotype.¶
In [4]:
data.shape
Out[4]:
We are checking how many good and bad predictions we have.¶
In [5]:
good_prediction = data[data['genotype'] == 1].shape[0]
bad_prediction = data[data['genotype'] == 0].shape[0]
In [6]:
print(str(round(good_prediction * 100 / (good_prediction + bad_prediction), 2)) + '% good SNP')
print(str(round(bad_prediction * 100 / (good_prediction + bad_prediction), 2)) + '% bad SNP')
Change genotype destingnation. We do that, becouse we want to search for 1 not for 0.¶
In [7]:
data['genotype'][data['genotype'] == 0.0] = 2.0
data['genotype'][data['genotype'] == 1.0] = 0.0
data['genotype'][data['genotype'] == 2.0] = 1.0
Becouse we read nucleotides as nucleotides triplets we merge those colmuns together as BEFORE and BEHIND¶
In [8]:
data['BEFORE'] = data['BEFORE1'] + data['BEFORE2'] + data['BEFORE3']
data['BEHIND'] = data['BEHIND1'] + data['BEHIND2'] + data['BEHIND3']
data = data.drop(['BEFORE1','BEFORE2','BEFORE3','BEHIND1','BEHIND2','BEHIND3'], axis=1)
In [9]:
data.profile_report()
Out[9]:
2. Preprocessing¶
2.1 Dealing with missing data.¶
Knowing that thare is missing data in variable DP2 and GQ we need to do something with that.¶
We have a few ways to deal with it :¶
- Delete those variables from our dataset
- Delete rows where is missing data
- We can do an imputation with mean, median, KNN, or even linear regression
- There is recommended way in deep learning to impute missing data with numer 0, or -9999, but when our variable is integer, or float it's not recommended then.
In our data we can see that in variable GQ there is 24.6% missing data, that's quite a lot. We can decide then to delete this variable from out dataset. In case of DP2 there is only 1% so we can make some imputation there. I decided to do it with mean, becouse KNN algorithm had some hardware problems.¶
In [10]:
data['DP2']=data['DP2'].fillna(data['DP2'].mean())
Deleting columns with missing data (GQ)¶
In [11]:
data = data.dropna(axis=1)
In [12]:
data = pd.get_dummies(data, prefix=['BEFORE', 'BEHIND'])
Now we can check how many variables we have now.¶
In [13]:
data.shape
Out[13]:
Quick look at our data now. We have nucleotides triplets from before and after variables, if a given nucleotide trio was in row then there appears integer 1.¶
In [14]:
data.head()
Out[14]:
2.2 Separation of data into X - features and Y - labels.¶
In [92]:
data_x = data.iloc[:, 1:]
data_y = data.iloc[:, :1]
2.3 Separation of data into test and train data. I decided to set test_size=0.1 becouse we have big data and I think that's enough.¶
In [93]:
x_train, x_test, y_train, y_test = train_test_split(data_x, data_y, stratify=data_y, test_size=0.1, random_state=10)
Check if there is good ratio good SNP to bad SNP:¶
In [94]:
y_train.genotype.value_counts()[0] / (y_train.genotype.value_counts()[1] + y_train.genotype.value_counts()[0])
Out[94]:
In [95]:
y_test.genotype.value_counts()[0] / (y_test.genotype.value_counts()[1] + y_test.genotype.value_counts()[0])
Out[95]:
2.4 Standardization with Z-score method ( subtract the average from the values and divide by standard deviation )¶
In [96]:
def standardize(column):
mean = x_train[column].mean()
std = x_train[column].std()
x_train.loc[:, column] = (x_train[column] - mean) / std
x_test.loc[:, column] = (x_test[column] - mean) / std
In [97]:
standardize('QUAL')
standardize('DP')
standardize('DP2')
standardize('CALL')
2.5 From traning data we can separate validation data.¶
In [98]:
x_train, x_val, y_train, y_val = train_test_split(x_train, y_train, stratify=y_train, test_size=0.1, random_state=11)
3. Building first, simple neural network.¶
We check what is the size/shape of first hidden layer¶
In [99]:
input_shape = x_train.shape[1]
x_train.shape
Out[99]:
In [23]:
def recall_m(y_true, y_pred):
true_positives = K.sum(K.round(K.clip(y_true * y_pred, 0, 1)))
possible_positives = K.sum(K.round(K.clip(y_true, 0, 1)))
recall = true_positives / (possible_positives + K.epsilon())
return recall
def precision_m(y_true, y_pred):
true_positives = K.sum(K.round(K.clip(y_true * y_pred, 0, 1)))
predicted_positives = K.sum(K.round(K.clip(y_pred, 0, 1)))
precision = true_positives / (predicted_positives + K.epsilon())
return precision
def f1_m(y_true, y_pred):
precision = precision_m(y_true, y_pred)
recall = recall_m(y_true, y_pred)
return 2*((precision*recall)/(precision+recall+K.epsilon()))
In [24]:
model = Sequential()
model.add(Dense(512, input_shape=(134,), activation='relu'))
model.add(Dense(256, activation='relu'))
model.add(Dense(128, activation='relu'))
model.add(Dense(64, activation='relu'))
model.add(Dense(32, activation='relu'))
model.add(Dense(16, activation='relu'))
model.add(Dense(8, activation='relu'))
model.add(Dense(1, activation='sigmoid'))
model.compile(optimizer='adam',
loss='binary_crossentropy',
metrics=[f1_m])
We can check params¶
In [25]:
model.summary()
This functions follow our NN and when F1 score in validation data will not imporve in 5 epochs in a row, then NN will stop.¶
In [26]:
callback_list = [
EarlyStopping(
monitor='f1_m',
patience=5,
mode='max'
),
ModelCheckpoint(
filepath='my_model.h5',
monitor='f1_m',
save_best_only=True,
mode='max'
)
]
3.3 Let's begin deep learning with our model.¶
In [27]:
history = model.fit(x_train, y_train,
batch_size=512,
epochs=200,
callbacks=callback_list,
validation_data=(x_val, y_val))
4. Results¶
4.1 Let's check first results¶
In [28]:
results = model.evaluate(x_test, y_test)
print('loss: ', results[0])
print('F1-score: ', results[1])
4.2 Confusion matrix¶
In [29]:
y_pred = model.predict(x_test)
y_pred_class = model.predict_classes(x_test)
cm = confusion_matrix(y_test,y_pred_class)
cr = classification_report(y_test,y_pred_class)
In [30]:
ax= plt.subplot()
sn.heatmap(cm, annot=True, fmt='g', cmap="Blues", ax = ax)
ax.set_xlabel('True labels');ax.set_ylabel('Predicted labels')
ax.set_title('Confusion Matrix')
ax.xaxis.set_ticklabels(['0', '1'])
ax.yaxis.set_ticklabels(['0', '1'])
Out[30]:
4.3 Sensitivity(Recall) and Precision¶
In [31]:
print(cr)
print('Precision:', str(round(cm[0][0] * 100 / (cm[0][0] + cm[1][0]), 2)) + '%')
print('Sensitivity:', str(round(cm[0][0] * 100 / (cm[0][0] + cm[0][1]),2)) + '%')
print('Specificity:', str(round((cm[1][1] *100 / (cm[1][1] + cm[0][1])),2)) + '%')
Sensitivity = True Positives / (True Positives + False Negatives)¶
Recall metric shows how many relevant samples are selected, which means how well our model can predict all the interested samples in our dataset.
Precision = True Positives / (True Positives + False Positives)¶
Precision metric tells us how many predicted samples are relevant i.e. our mistakes into classifying sample as a correct one if it’s not true.
Specificity = 1 - (False Positives / (False Positives + True Negatives))¶
Protability that for good SNP (0) classification will be correct (0)4.4 Balanced accuracy score¶
The balanced_accuracy_score function computes the balanced accuracy, which avoids inflated performance estimates on imbalanced datasets. It is the macro-average of recall scores per class or, equivalently, raw accuracy where each sample is weighted according to the inverse prevalence of its true class. Thus for balanced datasets, the score is equal to accuracy. In the binary case, balanced accuracy is equal to the arithmetic mean of sensitivity (true positive rate) and specificity (true negative rate), or the area under the ROC curve with binary predictions rather than scores. [https://scikit-learn.org/ ]¶
In [32]:
balanced_accuracy_score(y_test, y_pred_class)
Out[32]:
4.5 Matthews Correlation Coefficient¶
The Matthews correlation coefficient is used in machine learning as a measure of the quality of binary (two-class) classifications. It takes into account true and false positives and negatives and is generally regarded as a balanced measure which can be used even if the classes are of very different sizes. The MCC is in essence a correlation coefficient value between -1 and +1. A coefficient of +1 represents a perfect prediction, 0 an average random prediction and -1 an inverse prediction. The statistic is also known as the phi coefficient. In the binary case TP, TN, FP and FN are respectively the number of true positives, true negatives, false positives and false negatives, the MCC is defined as :¶
$ MCC = \frac{tp \times tn - fp \times fn}{\sqrt{(tp + fp)(tp + fn)(tn + fp)(tn + fn)}}. $
In [33]:
matthews_corrcoef(y_test, y_pred_class)
Out[33]:
5. Plots¶
5.1 Plot show difference between increase F1 score in train and test data. We can see how F1 score was changing threw epochs, that's really cool and coomon plot, but not necessary.¶
In [34]:
plt.plot(history.history['f1_m'])
plt.plot(history.history['val_f1_m'])
plt.title('model F1 score')
plt.ylabel('F1')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='upper left')
plt.show()
5.2 Plot - differance between decline in loss and validation loss.¶
In [35]:
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='upper left')
plt.show()
In [62]:
lw = 2
fpr, tpr, thresholds = roc_curve(y_test, y_pred)
roc_auc = auc(fpr, tpr)
def plot_roc_curve(fpr,tpr):
plt.plot(fpr,tpr,color='darkorange',
lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.axis([0,1,0,1])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.legend(loc="lower right")
plt.show()
plot_roc_curve (fpr,tpr)
6. Second Neural Network with weights.¶
The most important thing there is that I gave weights for our class. It's popular in inbalanced data, becouse we raise value of minority class.¶
- Large Weight: Assigned to examples from the minority class.
- Small Weight: Assigned to examples from the majority class.
I used this time dropout method to prevent overfitting, becouse we could see that early stopping isn't really good with our data.¶
In [100]:
model = Sequential()
model.add(Dense(512, input_shape=(134,), activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(256, activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(128, activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(64, activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(32, activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(16, activation='relu'))
model.add(Dense(1, activation='sigmoid'))
model.compile(optimizer='adam',
loss='binary_crossentropy',
metrics=[f1_m])
In [101]:
weights = {0:1, 1:5}
history = model.fit(x_train, y_train,
class_weight=weights,
batch_size=512,
epochs=100,
validation_data=(x_val, y_val))
7. Results¶
In [102]:
results = model.evaluate(x_test, y_test)
print('loss: ', results[0])
print('F1-score: ', results[1])
7.1 Confusion matrix¶
In [103]:
y_pred = model.predict(x_test)
y_pred_class = model.predict_classes(x_test)
cm = confusion_matrix(y_test,y_pred_class)
cr = classification_report(y_test,y_pred_class)
ax= plt.subplot()
sn.heatmap(cm, annot=True, fmt='g', cmap="Blues", ax = ax)
ax.set_xlabel('True labels');ax.set_ylabel('Predicted labels')
ax.set_title('Confusion Matrix')
ax.xaxis.set_ticklabels(['0', '1'])
ax.yaxis.set_ticklabels(['0', '1'])
Out[103]:
7.2 Other results, same as earlier¶
In [116]:
print(cr)
print('Precision:', str(round(cm[0][0] * 100 / (cm[0][0] + cm[1][0]), 2)) + '%')
print('Sensitivity(Recall):', str(round(cm[0][0] * 100 / (cm[0][0] + cm[0][1]),2)) + '%')
print('Specificity:', str(round((cm[1][1] *100 / (cm[1][1] + cm[0][1])),2)) + '%')
7.2.1 Balanced Accuracy score¶
In [105]:
balanced_accuracy_score(y_test, y_pred_class)
Out[105]:
7.2.2 MCC¶
In [106]:
matthews_corrcoef(y_test, y_pred_class)
Out[106]:
8. Plots¶
8.1 F1 score threw epochs¶
In [107]:
plt.plot(history.history['f1_m'])
plt.plot(history.history['val_f1_m'])
plt.title('model F1 score')
plt.ylabel('F1')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='upper left')
plt.show()
8.2 Loss threw epochs¶
In [108]:
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.title('model loss')
plt.ylabel('loss')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='upper left')
plt.show()
8.3 ROC plot¶
In [109]:
lw = 2
fpr, tpr, thresholds = roc_curve(y_test, y_pred)
roc_auc = auc(fpr, tpr)
plot_roc_curve (fpr,tpr)
9. Last Neural Network with StratifiedKFold.¶
In nutshell StratifedKFold take different part of our data K times for test data and every time there will be same percentage of genotype = 1 in every data. The data will be distributed equally.¶
I made there a ROC PLOT too, becouse that was easy way to create it and show all fold's curve.¶
In [110]:
skf = StratifiedKFold(n_splits=6, shuffle=True, random_state=3)
models_results = []
mean_tpr = 0.0
mean_fpr = np.linspace(0, 1, 100)
all_tpr = []
i = 1
for train_index, test_index in skf.split(data_x, data_y):
x_train, x_test = data_x.iloc[train_index], data_x.iloc[test_index]
y_train, y_test = data_y.iloc[train_index], data_y.iloc[test_index]
standardize('QUAL')
standardize('DP')
standardize('CALL')
standardize('DP2')
x_train, x_val, y_train, y_val = train_test_split(x_train, y_train, stratify=y_train, test_size=0.1, random_state=11)
model = Sequential()
model.add(Dense(512, input_shape=(134,), activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(256, activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(128, activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(64, activation='relu'))
model.add(Dropout(0.2))
model.add(Dense(32, activation='relu'))
model.add(Dense(16, activation='relu'))
model.add(Dense(8, activation='relu'))
model.add(Dense(1, activation='sigmoid'))
model.compile(optimizer='adam',
loss='binary_crossentropy',
metrics=[f1_m])
weights = {0:1, 1:7.5}
history = model.fit(x_train, y_train,
class_weight=weights,
batch_size=512,
epochs=25,
validation_data=(x_val, y_val))
loss_and_f1 = model.evaluate(x_test, y_test)
y_pred = model.predict(x_test)
y_pred_class = model.predict_classes(x_test)
cm = confusion_matrix(y_test,y_pred_class)
cr = classification_report(y_test,y_pred_class)
results = (loss_and_f1[0], loss_and_f1[1], cm, cr, y_pred_class, y_test)
models_results.append(results)
fpr, tpr, thresholds = roc_curve(y_test, y_pred)
mean_tpr += np.interp(mean_fpr, fpr, tpr)
mean_tpr[0] = 0.0
roc_auc = auc(fpr, tpr)
plt.plot(fpr, tpr, lw=1, label='ROC fold %d (area = %0.2f)' % (i, roc_auc))
i+=1
plt.plot([0, 1], [0, 1], '--', color=(0.6, 0.6, 0.6), label='Luck')
mean_tpr /= 6
mean_tpr[-1] = 1.0
mean_auc = auc(mean_fpr, mean_tpr)
plt.plot(mean_fpr, mean_tpr, 'k--',
label='Mean ROC (area = %0.2f)' % mean_auc, lw=2)
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic (ROC) with cross validation ')
plt.legend(loc="lower right")
plt.show()
10. Results¶
10.1 Results for every fold.¶
In [112]:
sr_acc = 0
sr_mcc = 0
print('Results for every fold: ' + '\n\n')
for i, j in enumerate(models_results):
print('-- Fold',i + 1,'--\n')
print('F1 score:',str(round(j[1], 4) * 100) + '%')
print('Loss:',round(j[0], 4))
print('\n' + 'Confusion matrix:' + '\n',j[2])
print('\n' + 'Precision:', str(round(j[2][0][0] * 100 / (j[2][0][0] + j[2][1][0]),2)) + '%')
print('Sensitivity(Recall):', str(round(j[2][0][0] * 100 / (j[2][0][0] + j[2][0][1]), 2)) + '%')
ba_acc = balanced_accuracy_score(j[5], j[4])
print('\n' + 'Balanced Accuracy Score:(=AUC)', ba_acc )
sr_acc += ba_acc
mcc = matthews_corrcoef(j[5], j[4])
print('\n' + 'MCC:', mcc)
sr_mcc += mcc
print('\n',j[3])
print('\n\n')
10.2 Averge results¶
10.2.1 Averge confusion matrix.¶
In [113]:
sm = models_results[0][2] + models_results[1][2] + models_results[2][2] + models_results[3][2] + models_results[4][2] + models_results[5][2]
av_cm = sm//6
In [114]:
ax= plt.subplot()
sn.heatmap(av_cm, annot=True, fmt='g', cmap="Blues", ax = ax)
ax.set_xlabel('True labels');ax.set_ylabel('Predicted labels')
ax.set_title('Confusion Matrix')
ax.xaxis.set_ticklabels(['0', '1'])
ax.yaxis.set_ticklabels(['0', '1'])
Out[114]:
10.2.2 Averge results for metrics.¶
In [121]:
print('Averge results:')
average = str(round((models_results[0][1] + models_results[1][1] + models_results[2][1] + models_results[3][1] +
models_results[4][1] + models_results[5][1]) / 6, 4)*100) + '%'
print('\n' + 'F1 score:', average)
print('Precision:', str(round(av_cm[0][0] * 100 / (av_cm[0][0] + av_cm[1][0]), 2)) + '%')
print('Sensitivity(Recall):', str(round(av_cm[0][0] * 100 / (av_cm[0][0] + av_cm[0][1]),2)) + '%')
print('Mean balanced accuracy score(=AUC):', str(round(sr_acc / 6, 4)))
print('Mean MCC:', str(round(sr_mcc / 6, 4)))
11. Summary¶
| Recall | Precision | balanced accuracy score | MCC | F1-SCORE | |
|---|---|---|---|---|---|
| Frist NN | 98.66% | 97.89% | 0.68 | 0.40 | 23% |
| Second NN | 99.08% | 97.95% | 0.69 | 0.46 | 25% |
| Last NN (averge results) | 97.52% | 98.21% | 0.7235 | 0.41 | 13.93% |
12. Conclusions¶
- Data I was working on is not easy to deal with. There should be a lot more work to do with this, and variables do not describe good and bad genotype well.
- Best result gave second Neural Network with weigths and dropout method.
- In first model occur overfitting, so results wasn't really good.
- Low F1 score in third NN is due to high weigths I used.
- StratifiedKfold could be a good option, but it takes a lot of time to calculate.
References¶
- Training deep neural networks on imbalanced data sets, Authors: Shoujin Wang, Wei Liu, Jia Wu, Longbing Cao, Qinxue Meng, Paul J. Kennedy (2016r.)
- https://www.kaggle.com/parulpandey/10-simple-hacks-to-speed-up-your-data-analysis
- https://www.scipy.org/
- https://scikit-learn.org/stable/
- https://machinelearningmastery.com/cost-sensitive-neural-network-for-imbalanced-classification/